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ABSTRACT 


In  processing  networks,  ordinary  network  constraints  are  supplemented 
by  proportional  flow  restrictions  on  arcs  entering  or  leaving  some 
nodes.  This  paper  describes  a  new  primal  partitioning  algorithm  for 
solving  pure  processing  networks  using  a  working  basis  of  variable 
dimension.  In  testing  against  MPSX/370  on  a  class  of  randomly  generated 
problems,  a  FORTRAN  implementation  of  this  algorithm  was  found  to  be 
an  order  of  magnitude  faster.  Besides  indicating  the  use  of  our  methods 
in  stand  alone  fashion,  the  computational  results  also  demonstrate 
the  desirability  of  using  these  methods  as  a  high-level  module  in  a 
mathematical  programming  system. 
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Pure  processing  networks  are  minimum  cost  flow  problems  in  which 
proportional  flow  restrictions  are  permitted  on  the  arcs  entering  or 
leaving  a  node  Applications  are  widespread  with  the  proportional  flow 
restrictions  governing  such  things  as  the  size  of  loan  payments  in  cash 
flow  models  and  the  relation  between  raw  materials  and  finished 
products  in  assembly  models  A  survey  of  applications  is  included  in 
the  work  of  Koene  [16]  The  proportional  flow  restrictions  can  be 
modeled  either  as  nonnetwork  rows  or  as  nonnetwork  columns  in  a  linear 
programming  ( LP )  formulation  Our  approach  uses  nonnetwork  columns 
since  they  lead  to  an  LP  basis  with  fewer  rows  In  [  1 6  ]  ,  Koene  shows 
that  any  LP  problem  can  be  readily  transformed  to  a  pure  processing 
network  problem  at  the  expense  of  enlarging  the  problem  size.  However, 
the  primal  partitioning  methods  of  this  paper  will  only  be  more 
efficient  than  standard  LP  methods  when  the  basic  nonnetwork  columns 
form  a  small  fraction  of  all  basic  columns  The  allowable  size  of  this 
fraction  is  vet  to  be  determined 

The  success  of  primal  simplex  solution  procedures  for  solving  pure 
networks  is  well  known  These  procedures  depend  on  special  data 
structures  popularized  by  Glover  and  klingman  and  t  he i r  co-workers  more 
than  a  decade  ago  [7],  [8],  [  1 3  ]  A  detailed  description  of  network 
methods  is  given  by  Bradley  et  al  [2]  Standard  primal  simplex  LP 
codes,  however,  use  data  structures  for  exploiting  sparsity  in  the 
basis  matrix  [20], [22]  As  shown  in  [8],  [9],  [25],  specialized  net¬ 
work  codes  have  achieved  an  advantage  in  solution  speed  up  to  two 
orders  of  magnitude  over  standard  LP  codes 


-2- 


3 


>8 


8 
* > 


Since  LP  problems  often  have  a  large  network  component,  ways  to 
exploit  this  component  based  on  specialized  network  methods  have  been 
sought.  This  creates  a  need  to  reconcile  the  two  data  structure  types 
Two  levels  of  detail  for  combining  these  data  structures  have  been 
used.  Glover  et  al  ,  [12]  describe  a  high-level  approach  in  which  a 

PL/1  main  program  calls  both  a  FORTRAN  network  code  and  MPSX'370  [14} 
modules  in  a  complex  solution  procedure  for  a  large  nonlinear  mixed 
integer  program.  McBride  [id]  and  Glover  and  Klingman  |lO].[ll] 
describe  solution  methods  for  embedded  network  problems  in  which  basis 
partitioning  allows  both  network  and  sparse  matrix  data  structures 
which  are  tightly  integrated  in  maintaining  the  basis.  The  specialized 
FORTRAN  code  of  [18]  ran  about  five  times  faster  than  MINOS  [21  ]  In 
[  1  I  ]  .  three  large  problems  which  included  both  nonnetwork  rows  and 
columns  were  solved  using  a  FORTRAN  embedded  network  code  and  MPSX./370 
The  nonnetwork  rows  make  these  problems  more  difficult  than  the  test 
problems  of  the  present  paper,  and  MPSX  '370  solved  them  about  4"  faster 
than  the  specialized  code  These  efforts  provide  a  start  in  delimiting 
the  class  of  problems  which  benefit  from  tight  integration  of  the  two 
data  structures 

Tomlin  and  Welch  [25]  describe  a  ma t  tiema t 1 < a  1  programming  system 
written  in  assembly  language  which  contains  two  optimizers,  one  based 
on  network  data  structures  and  one  based  on  sparse  matrix  data 
structures  Some  modifications  of  network  methods  were  made  in  order 
to  accommodate  the  network  optimizer  into  the  system  This  is  a 
high-level  approach  since  problems  which  are  partly  network  do  not 
benefit  from  the  specialized  network  data  structures  However,  the  two 
optimizers  do  have  common  I/O  and  start  routines.  For  embedded  network 
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problems  with  nonnetwcrk  rows,  a  two  stage  starting  procedure  is 
described  in  which  the  network  portion  of  the  problem  is  first  solved 
using  the  network  optimizer  This  solution  is  then  used  to  provide  a 
partial  basis  for  the  regular  LP  optimizer  Presumably  this  approach 
is  superior  to  using  only  the  regular  LP  optimizer  and  its  start 
routine,  however,  comparative  solution  times  are  not  provided  in  [25] 

A  similar  approach  was  used  in  [il]  where  a  FORTRAN  network  optimizer 
solved  Lagrangian  relaxations  of  the  original  problem  iteratively  to 
provide  an  initial  partial  basis  for  MPSX,  370  This  start  procedure 
resulted  in  a  total  solution  time  which  was  21  times  smaller  than  that 
achieved  when  the  MPSX/370  CRASH  start  procedure  was  used 

Chen  and  Engquist  [5]  described  a  primal  simplex  variant  which  is 
the  precursor  of  the  algorithm  of  this  paper  One  feature  of  the 
algorithm  of  [5]  is  that  all  nonnetwork  columns  are  always  basic.  This 
results  in  a  working  basis  whose  dimension  is  equal  to  the  number  of 
nonnetwork  columns  When  a  FORTRAN  implementation  of  the  previous 

algorithm  was  tested  against  MINOS,  it  was  found  to  be  an  order  of 
magnitude  faster  on  problems  with  up  to  2UU  nonnetwork  columns  The 
question  remained  as  to  whether  a  similar  improvement  in  efficiency 
could  be  achieved  against  an  assembly  language  code  such  as  MPSX/370 
In  the  present  paper,  we  answer  this  question  in  the  affirmative  We 
describe  improvements  to  both  the  algorithm  and  implementation  of  [5] 
In  particular,  the  working  basis  for  the  current  method  has  a  dimension 
which  varies  with  the  number  of  basic  nonnetwork  columns  The  approach 
we  use  involves  a  fairly  tight  integration  of  network  and  sparse  matrix 
data  structures.  The  design  of  our  processing  network  code,  PROCNET, 
was  influenced  by  Marsten  [17]  in  that  it  consists  of  a  library  of 


subroutines  which  communicate  through  parameter  lists.  Like  Marsten's 
XMP  code,  PROCNET  utilizes  the  LA05  subroutines  of  Reid  [23],  [24]  In 


fact,  our  extension  of  the  LA05  subroutines  may  be  beneficial  to  XMP 
users.  This  is  discussed  in  Section  4  Although  PROCNET  is  already 

highly  efficient,  further  speedup  is  no  doubt  possible  through  a 
tighter  integration  of  network  and  sparse  matrix  (LA05)  data 
structures  It  would  also  be  of  interest  in  future  work  to  use  PROCNET 
as  a  module  in  a  high-level  approach.  For  example,  in  large  embedded 
network  problems  having  both  nonnetwork  constraints  and  nonnetwork 
variables.  the  nonnetwork  constraints  could  be  initially  relaxed 
PROCNET  could  then  be  used  to  solve  the  remaining  problem  and  thus  pro¬ 
vide  a  partial  starting  basis  for  an  LP  optimizer. 

2  Background  on  Processing  Networks 

The  pure  processing  network  problem  (problem  PPN)  is 


mimmi  ze 

CNXN  + 

Cp>p 

(2 

1) 

subject  to 

aNxN  + 

A,,Xp  -  b 

(2 

21 

o  $  XN 

$  hN 

(2 

3) 

{)  ^  Xp 

(2 

4  1 

The  mxn  matrix  is  the  node  arc  incidence  matrix  for  a  pure  net¬ 
work  N.  while  the  mxp  matrix  Ap,  contains  the  nonnetwork  columns  These 
nonnetwork  columns  are  also  called  processing  columns.  Vector  b 
contains  the  supply  values,  while  c^  and  Cp  contain  unit  costs  for  the 
vectors  of  decision  variables  x^  and  Xp  The  capacities  are  the 
components  of  the  simple  upper  bound  vectors  h^j  and  hp  It  is  assumed, 
without  loss  of  generality,  that  a  slack  arc  and  artificial  arcs  with 


Big-M  costs  have  been  introduced  into  the  network  N  so  that  it  is 
connected  and  the  matrix  has  rank  m.  It  is  also  assumed  that  each 
row  of  [A^.Ap]  contains  at  most  one  non-zero  component  from  the  columns 
of  Ap.  The  latter  assumption  is  made  to  simplify  the  notation,  and  it 
does  not  restrict  the  application  of  our  methods.  Each  column  of  Ap  is 
of  the  form 

a  in  row  v 

-«vw(z)  ln  row  w*z>'  2=1,2 . t  (2.5) 

0  e 1 sewhere 

Further,  av  =1,0*;  a  «;  1.  and 

t 

I  avw(  z)  “  1  (26) 

z  =  l 

must  hold 

Corresponding  to  each  column  of  Ap  is  a  column  of  A^  called  an  alloca¬ 
tion  column  The  following  network  diagram  is  associated  with  this 


struct  ure 


Figure  1  A  Splitting  Node 

In  Figure  1,  the  square  is  termed  a  splitting  node,  while  the  circles 
are  ordinary  network  nodes  Arc  (u.v)  is  called  an  allocation  arc  and 
its  column  in  is  the  allocation  column  The  arcs  leaving  the  square 
node  in  Figure  1  are  termed  processing  arcs,  and  they  are  represented 
by  a  column  of  the  form  (2  f> )  Node  v  and  nodes  w(  1  )  .  ,w(t)  in 
Figure  l  are  called  processing  nodes  Conservation  of  flow  and  the 
conditions  on  the  avw  imply  that  a  flow  x  entering  the  splitting  node 


in  Figure  1  generates  a 

f low  ftvw(z)x  alone 

processing 

arc  ( v ,w( z )  ) 

To 

be  consistent  with  (2 

5)  we  assume  that 

ord l nary 

ne  twork 

arcs 

are 

represented  by  a  column  of  A^  which 

contains 

a  1  in 

the 

row 

corresponding  to  the 

ta i 1  node  of  the 

arc  and 

a  -1  in 

the 

row 

corresponding  to  the  head  node  In  Figure  1,  if  arc  (u.v)  corresponds 
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to  column  r  of  A^j  and  the  corresponding  processing  column  (2.5)  is 
column  s  of  Ap.  then  it  is  assumed  that  the  capacities  [h^]r  and  [hp]s 
are  equal . 

In  [16],  the  definition  of  pure  processing  networks  includes  the 
structure  formed  when  the  direction  of  the  arcs  in  Figure  1  is 

reversed  For  problems  with  finite  capacities,  by  complementing  flows 
with  respect  to  these  capacities  and  adjusting  supply  values 
appropriately,  this  structure  can  be  reduced  to  the  one  shown  in 
Figure  1  Thus,  there  is  no  loss  ot  generality  in  restricting 

attention  to  the  structure  of  Figure  1 

In  order  to  exploit  the  network  structure  of  1’PN ,  it  is  necessary 
to  see  how  this  structure  carries  over  to  a  basis  The  first 

observation  to  be  made  is  that  the  slack  arc  column  must  be  present  in 
any  basis  matrix,  otherwise  the  rows  of  the  matrix  would  sum  to  zero 
Next,  we  note  that  when  the  processing  columns  and  the  slack  column  are 
removed  from  a  basis  matrix  containing  r  processing  columns,  the 
remaining  columns  correspond  to  the  arcs  in  r+1  trees  This  follows  by 
a  simple  counting  argument  These  r+1  trees  are  called  basis  trees 

The  basis  tree  which  is  i nc  ident  to  the  problem's  slack  arc.  when  taken 
together  with  that  arc  .  forms  the  basis  quasi-tree  For  the  remaining 
r  trees,  root  nodes  are  chosen  arbitrarily  and  the  resulting  r  rooted 
trees  are  called  the  rooted  basis  trees  We  let  a  basis  matrix  B 

containing  r  processing  columns  be  partitioned  as 

B  =  [B, .Bg]  (2  7) 


'■  *’  .  «**/■  l  .  A.  V* 

» *  j.  ►  j»  * jkY*  r m  \m  v .  •- w  .  • 


<** ^ W.  **-  *  . 
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where  Bj  contains  network  columns  and  Bg  contains  the  processing 


columns  . 


If  the  last  r  rows  of  B  correspond  to  roots  of  rooted  basis  trees, 


T  C 

and  = 

D  F 


(2  8) 


results  where  T  and  C  have  m-r  rows  and  D  and  F  have  r  rows  The 
working  basis  corresponding  to  basis  B  is  defined  to  be  the  matrix  Q 

where 


Q  =  F  DT  *C 


(2.9) 


It  can  be  shown  that  Q  is  nonsingular,  see,  e  g  ,  [5],  [ 1 5 ] 

The  matrix  T  in  (2.8)  corresponds  to  a  collection  of  quasi-trees 

By  judicious  use  of  matrices  T  and  Q,  updated  entering  columns  and  dual 

variables  for  the  primal  simplex  algorithm  can  be  computed  without  the 

need  for  maintaining  a  factorization  of  the  entire  basis  matrix  B  If 

a  =  [  1  ]  is  the  entering  column  and  it  is  partitioned  to  be 

a2 

compatible  with  (2.8),  then  a  straightforward  calculation  shows  that 

y 

the  updated  entering  column  [  1  ]  is  obtained  by  solving 

*2 


Ov^  =  an  -  DT  a, 


(2.10) 


Similarly,  we  partition  the  basic  costs  Cg^Cj.Cg]  and  the  dual 
variables  rr=  f  rr  j  ,  1  so  that  they  are  compatible  with  (2  8).  The  dual 
variables  are  computed  by  solving 

n2Q  =  c2  -  c  j T_ 1 C  (2  12) 

n ,  =  CjT-1  -  n.-,DT~ 1  (213) 


Suppose  that  column  j  of  B?  is  the  processing  column  with 
splitting  node  v  p  is  defined  to  be  the  set  of  processing  nodes  in 
rooted  basis  tree  i  which  correspond  to  nonzero  values  m  column  j  of 
Bp  The  following  theorem  was  proved  in  [5) 

Theorem  For  a  basis  H.  the  elements  q  (  of  the  working  basis  Q 
sat  i s  f  v 


where  the  sum  in  (2  14)  is  defined  to  be  zero  in  case  P  is  empty 

If  tin  t  li  t  e  l  i  lig  t  w  I  in  in  i  i  ii  i  r  es  pul  ids  to  till  Ml  i  with  both  end  nodes 
iri  a  common  basis  tree  it  follows  from  (  3  10)  and  (2  11)  that  the  only 
flow  changes  occur  on  the  cycle  in  the  tree  formed  bv  the  entering  arc 
Furthermore.  the  theorem  implies  that  no  working  basis  update  is 
required  in  this  case  This  type  of  pivot  is  termed  a  pure  network 


pivot  wti  i  1 e  all  other  pivots  are  termed  processing  network  pivots 


3  Primal  Simplex  Variant 


The  fundamental  observation  in  the  development  of  the  primal 
simplex  variant  which  we  use  relates  to  Figure  1  The  flow  on  the  al¬ 
location  arc  (u.v)  must  equal  thp  value  of  the  variable  for  the 
associated  processing  column  Thus,  if  the  processing  column  is  to  be 
the  leaving  variable  the  n i  I m  a  I  i on  an  can  leave  the  basis  instead 
Note  that  the  all  o(  a  I  i  on  si  i  must  tie  basic  in  this  situation,  since 
otherwise  the  pivot  would  lead  to  the  impossible  situation  in  which 
both  the  allocation  arc  and  the  processing  column  are  nonbasic. 

Before  stating  the  simplex  algorithm,  we  outline  the  situations 
which  are  to  be  considered  in  updating  the  basis  trees  and  the  working 
basis  Q  during  the  basis  exchange  step  Before  the  basis  exchange  is 
executed,  we  assume  that  is  the  basis  quasi-tree  and  t  ,  1=1,2,  ,r 
are  the  rooted  basis  trees  Those  basis  trees  which  have  been  changed 
during  the  exchange  step  will  be  designated  by  means  of  an  asterisk. 
If  a  change  to  r(,  i^o,  results  in  a  change  to  one  of  the  sets  P  j  ^  in 
(2  14).  then  row  i  of  Q  must  be  updated 

Several  cases  occur  in  the  basis  exchange  step  of  the  simplex 
algorithm  However.  the  variant  we  describe  allows  us  to  restrict 
attention  to  the  two  following  cases 

i)  The  entering  column  is  a  processing  column  and  the  leaving  column 
is  a  network  column  (arc)  If  the  leaving  arc  is  in  r  ,  then  row 
l  of  Q  wi  1  1  be  updated  (unless  1=0)  and  an  additional  row  and 
column  will  be  adjoined  to  Q. 


(ii)  Both  the  entering  and  leaving  columns  are  network  columns  (arcs) 
If  the  entering  arc  is  incident  to  and  r  ,  then  these  two  trees 
are  joined  to  form  r*  .  If  the  leaving  arc  is  contained  in  . 
then  splits  into  two  trees  upon  its  removal  One  of  these 
becomes  r*  and  the  other  becomes  t*  If  1  ,  j  ,  and  k  are  nonzero 
mid  distinct,  then  three  rows  of  Q  will  be  updated  Otherwise, 
special  cases  occur  in  which  at  most  two  rows  of  Q  are  updated 
One  of  these  special  cases  is  the  pure  network  pivot  for  which  no 
updating  of  rows  of  Q  is  necessary 

We  remark  that  cases  in  which  a  processing  column  leaves  the  basis  are 

not  considered  here,  since  the  allocation  arc  can  be  chosen  to  leave 

i ns  t  ead 

The  lias  is  trees  can  be  visualized  as  hanging  downward  from  their 
roots  The  node  incident  to  the  slack  arc  in  the  basis  quasi-tree  is 
taken  as  the  root  there  If  two  basis  trees  and  are  joined  by  an 
entering  arc.  the  resulting  will  retain  the  root  of  r(  while  will 
banc  below  r  in  r*  Also,  when  a  leaving  are  is  deleted  from  a  basis 

lice  r  j.  .  an  upper  tree  i  j.  j  which  contains  the  root  of  11  lower 

t ree  Tg  ,  are  formed 

A  starting  F’f'N  basis  can  be  obtained  by  first  setting  the 
variables  Xp  in  (2  l)-(2  •!  1  to  upper  or  lower  bounds  according  to  some 
heuristic  procedure  These  variables  are  set  nonbasic  in  the  PPN 
basis  Those  >p  variables  at  upper  bounds  induce  supplies  in  network  N 
in  addition  to  those  represented  by  b  in  (22)  The  resulting  network 
problem  is  solved  to  optimality  to  give  the  initial  collection  of  basis 
trees  which,  in  this  case,  consists  of  a  single  quasi-tree  Of  course, 
this  quasi-tree  may  contain  artificial  arcs  with  positive  flow  An 


extension  of  this  starting  procedure  has  been  implemented  as  described 
in  Sect  ion  4 . 

We  introduce  the  vector  X  to  represent  certain  quantities  which 
may  be  thought  of  as  pseudo  node  potentials. 

X  =  CjT-1  (31) 

It  will  be  useful  to  extend  X  by  defining  X  =0  for  root  nodes  j  of 
rooted  basis  trees  For  convenience.  this  extension  will  also  be 
denoted  as  X 

I’rimal  Simplex  A  Igor  it  tun  for  PPN 

U  Obtain  an  initial  basis  Set  up  the  initial  basis  trees  and 
working  hasis  Compute  initial  dual  variables  and  basic  solution 
1  Price  nonbasic  processing  columns  and  arcs.  If  an  entering  proces¬ 
sing  column  is  found,  go  to  Step  3  If  an  entering  arc  is  found, 
go  to  Step  2.  I  f  no  entering  processing  column  or  arc  exists, 
check  for  basic  artificial  arcs  with  positive  flow  If  basic 
artificials  with  positive  flow  exist,  s t op  with  an  infeasible  prob¬ 
lem.  otherwise,  stop  with  an  optimal  solution 

If  both  end  nodes  of  entering  art  e  are  not  in  a  common  basis  tree 
t,  go  to  Step  3  Otherwise,  restrict  the  ratio  test  and  flow 
update  to  the  arcs  on  the  cycle  formed  in  r  by  e  Update  X  on  the 
tree  hanging  below  e  after  the  leaving  arc  is  removed  Go  to 


Compute  and  yg  using  (2.10)  and  (2.11). 

Perform  the  ratio  test  restricted  to  arcs.  Update  basic  solution 
values 

Update  basis  trees  and  working  basis  (basis  exchange  step)  If  an 
arc  is  entering  the  basis  go  to  (11) 

(i)  A  processing  column  enters  the  basis,  and  the  leaving  arc  is 

in  The  leaving  art  is  removed  to  form  an  upper  tree  Tj_j 

and  u  lower  tree  'j.p  Tree  becomes  rr+]'  ^  i  s  updated  on 

r  j.  and  row  r+1  of  Q  is  created  using  (2.14)  Column  r+1 
of  y  corresponding  to  the  entering  column  is  also  created 
using  (2  14)  Tree  ;  j.  t  becomes  r  £  and  row  k  of  Q  is  updated 
(unless  k~0)  Go  to  Step  f> 

(ii)  All  art  e  enters  the  basis  (The  details  follow  for  this 

step  when  e  is  incident  to  t(  and  Tjt  the  leaving  arc  is  in 

r.  .  and  i,  j.  k  are  nonzero  and  distinct  The  remaining 
Iv 

cases  involve  at  most  two  rows  of  Q  and  the  details  are 

omitted  )  First,  hangs  below  7(  via  arc  e  to  form  tJ  and  * 

is  updated  on  How  i  of  U  is  updated  to  form  Q’  Next  . 

the  leaving  art  is  removed  from  to  form  an  upper  tree  j 

and  a  lower  tree  r.  .,  The  lower  tree  becomes  t*  and  A  is 

K ..  j 

updated  oil  r  j  How  j  of  Q’  is  updated  to  form  Q*'.  Finally. 
t^j  becomes  T|*  and  row  k  of  Q' *  is  updated  to  form  Q*** 

Update  rig  using  (2  12).  Compute  rij  using 

n  j  =  A  —  rigDT-1  (3.2) 

where  A  has  been  previously  updated  Go  to  Step  1 


In  the  above  algorithm,  processing  columns  are  allowed  to  enter 
the  basis,  but  not  to  leave.  However,  when  the  basis  is  reinverted, 
eligible  processing  columns  are  removed  from  the  basis  via  a  series  of 
degenerate  pivots  More  details  on  this  procedure  are  given  in 
Section  4. 

Updating  of  X  on  a  tree  which  is  rehung  is  done  by  adding  a 
certain  constant  to  these  X  values  in  the  same  way  as  node  potentials 
are  updated  in  the  pure  network  cose. 

4  Implementation 

A  FORTRAN  implementation,  PROCNET ,  of  the  primal  simplex  algorithm 
of  Section  3  was  created  This  version  of  PROCNET  extends  and  enhances 
a  previous  version  which  is  described  in  [5]  Problem  data  storage  in 
PROCNET  is  accomplished  by  means  of  arrays  for  arc  costs,  capacities, 
and  head  nodes  Also,  arrays  containing  the  nonzero  values  in  proces¬ 
sing  columns  and  the  positions  of  these  values  are  used  The  costs  of 
processing  columns,  components  of  Cp  in  (2.1),  are  assumed  to  be  zero, 
since  such  costs  can  be  placed  on  the  allocation  arcs  instead 

The  basis  trees  are  incorporated  into  a  single.  larger  tree 
following  [10],  [ll]  This  tree  is  referred  to  as  the  master  basis 
tree  and  its  root  is  called  the  master  root  The  roots  of  all  basis 
trees  are  connected  to  the  master  root  by  artificial  arcs,  and  the 
slack  arc  of  the  basis  quasi-tree  is  disregarded  since  it  plays  no  role 
in  the  implementation  For  maintaining  the  master  basis  tree,  the 
predecessor,  depth,  thread  and  reverse  thread  functions  [2],  [13],  [15] 


are  emp 1 oyed 


-15- 


Since  the  updates  to  Q  in  Step  5  of  the  primal  simplex  algorithm 

involve  more  changes  to  rows  than  to  columns,  we  have  elected  to  main- 

T 

tain  Q  by  applying  column  replacement  techniques  to  its  transpose  Q 

Three  LA05  subroutines  were  written  in  FORTRAN  by  Reid,  and  they 
are  described  in  [23],  [24).  These  subroutines  implement  a  sparse 

variant  of  the  Bar t e  1  s-Go  1  ub  algontlim  |l]  In  order  to  utilize  these 
subroutines  for  maintaining  the  working  basis  for  PROCNET ,  we  needed  to 
extend  them  by  providing  a  means  for  adjoining  additional  rows  and 
columns.  The  two  subroutines  we  created  for  this  purpose  are  described 
in  this  section  and  we  note  that  they  may  be  useful  in  situations  other 
than  maintaining  a  working  basis  For  example,  Marsten's  XMP  linear 
programming  code  [17]  uses  the  original  LA05  subroutines  for  main¬ 
taining  the  LP  basis  Our  additional  subroutines  can  be  used  to  pro¬ 
vide  a  restart  capability  for  XMP  when  one  or  more  rows  are  adjoined  to 
an  LP  problem 

We  proceed  with  an  explanation  of  the  functions  of  the  three 

T 

original  I.A05  subroutines  as  applied  to  the  rxr  matrix  Q  in  PROCNFT 
The  LA05A  subroutine  produces  a  factorization 


Q1  =  LH 


(4  I  ) 


The  lower  triangular  factor  L  is  maintained  as  a  sequence  of  eta 
vectors  The  upper  triangular  factor  l1  is  stored  as  a  sparse  matrix 
with  the  nonzeros  in  the  rows  held  in  packed  form,  while  only  the 
positions  for  nonzeros  in  columns  are  kept  Additional  information 
which  is  maintained  includes  the  pivot  order  and  the  number  of  nonzeros 


per  row  or  column 


The  LA05B  subroutine  solves  sets  of  equations 


QT5  =  E 


(4.2) 


Q  x  =  E 


(4  3) 


where  x  and  E  are  r-d mens  1 onal  vectors  This  solution  is  carried  out 
with  the  use  of  (4.1). 

The  third  subroutine,  LA05C,  revises  the  factorization  (41)  when 
T 

one  of  the  columns  of  Q  is  changed 

T 

In  order  to  accommodate  additional  rows  and  columns  for'  Q  in 

— T 

(41),  we  embed  this  matrix  in  a  larger  matrix  Q  where 


1  0 


0  Q1 


(4.4) 


and  I  is  on  identity  matrix  of  dimension  s.  When  an  additional  row  d  of 
*  of  length  Ml  is  treated,  it  is  embedded  in  a  row  of  length  s  +  r  con¬ 
taining  s-1  initial  zeros  as 


[  0  ,  .  0  ,  d  ] 


(4.5) 


and  the  row  of  (4.5)  is  inserted  as  row  s  of  Q1  in  (44)  Likewise,  an 

T 

additional  column  of  Q  is  supplied  with  s-1  initial  zeros  and  inserted 
as  column  s  in  (44) 


Crii 


A  new  subroutine  LA05SS  was  created  which  takes  the  original  fac- 
T 

torization  of  Q  from  (4.1)  as  provided  by  LA05A  and  adjusts  the 
information  for  storing  L  and  U  to  produce  the  factorization 


where 


(46) 

(4  7) 


(4  0) 


Essent l a  1  1  y 

wha  t 

1  s 

done  is  to 
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A  new  subroutine  LA05TT  carries  out  the  task  of  inserting  a  new 
row  (45)  when  this  is  required  We  note  that  changes  are  needed  for 

the  data  maintaining  0  but  eta  vectors  corresponding  to  L  remain 

_t 

correct  For  insertion  of  n  new  <olumn  to  i)  I.A0f>C  is  used 

PROCNET  obtains  an  initial  basis  for  PPN  by  means  of  a  heuristic 
based  on  [3],  (Oj  In  order  to  apply  the  heuristic,  the  arc  data  for 
each  processing  arc  are  generated  The  resulting  pure  network  with 
proportional  flow  restrictions  relaxed  is  solved  first.  Next,  the  flow 
values  of  the  relaxed  solution  on  the  allocation  arcs  are  used  to 
creeate  a  new  pure  network  problem  with  nonzero  lower  bounds  on  the 
processing  arcs  If  the  flow  value  on  the  allocation  arc  (u.v)  of 
Figure  1  is  x,  then  the  lower  bound  on  arc  (v,w(z))  is  set  to 
0  75«vw(z)*  This  second  network  problem  is  solved  and  the  flows  on 


allocation  arcs  are  saved  If  such  flows  are  at  the  original  problem 
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bounds  ,  the  corresponding  processing  columns  are  made  nonbasic  while 
the  allocation  arc  is  basic  If  an  allocation  arc  has  a  flow  between 
the  original  problem  bounds,  a  parallel  allocation  arc  is  created  with 
a  capacity  equal  to  this  flow  value  The  parallel  arc  is  given  the 
same  cost  as  the  original  allocation  arc,  while  the  original  allocation 
arc  is  replaced  bv  a  modified  allocation  arc  whose  capacity  equals  the 
original  capacity  less  the  flow  value  In  the  initial  PPN  basis,  the 
modified  allocation  arc  is  nonbasic  at  zero,  the  parallel  arc  is 
nonbasic  at  capacity  and  the  corresponding  processing  column  is  basic 
The  initial  working  basis  Q  is  an  r*r  identity  matrix  where  r  is  the 
number  of  parallel  arcs  created  The  initial  flows  through  allocation 
arcs  and  their  parallel  arcs  induce  supply  values  on  the  remainder  of 
the  network  Tin s  remaining  network  problem  is  solved  to  optimality 
and  the  resulting  optimal  tree  becomes  the  initial  PPN  basis 
quas l -t  r ee 

We  note  that  an  improvement  to  the  way  allocation  arcs  and  their 
parallel  arcs  are  handled  has  resulted  in  about  a  10"  decrease  in  the 
number  of  PPN  iterations  over  the  previous  version  of  PROCNET  |5] 
Tin  s  improvement  is  accomplished  by  reinstating  the  original  allocation 
arc  and  eliminating  the  modified  allocation  arc  and  its  parallel  arc 
once  one  of  the  latter  arcs  enters  the  basis  In  the  previous  version 
of  the  code  both  the  modified  allocation  arc  and  its  parallel  arc  were 
maintained  for  all  PPN  pivots,  and  this  restricted  the  amount  of  flow 
change  which  could  be  achieved  on  a  single  pivot  involving  these  arcs. 


PROCNET  uses  two  conditions  to  trigger  a  basis  reinversion 


The 


first  of  these  conditions  is  a  total  of  40  column  and  row/column 
updates  of  Q  The  other  condition  involves  the  dimension  s  of  I  in 
(4.8).  After  s  processing  columns  have  entered  the  basis,  the  next 
entering  processing  column  causes  a  basis  reinversion  during  which  a 
new  identity  I  is  created  PROCNET  currently  uses  a  value  of  10  for  s 
At  the  time  of  basis  reinversion.  PROCNET  searches  for  basic  processing 
columns  having  corresponding  allocation  arcs  (see  Figure  1)  which  are 
nonbasic.  Before  this  basis  reinversion  takes  place,  a  series  of 

degenerate  pivots  is  executed  in  which  such  processing  columns  are  made 
nonbasic  while  their  allocation  arcs  become  basic 

Pricing  for  PROCNET  is  handled  by  means  of  two  candidate  lists.  LI 
and  1,2  LI  is  used  for  pure  network  pivots  while  processing  network 
pivots  arising  from  either  entering  processing  columns  or  arcs  are 
placed  on  L2  In  order  to  identify  arcs  which  correspond  to  pure  net¬ 
work  pivots,  basis  trees  are  numbered  A  node  length  array,  treenum. 
assigns  to  each  node  of  a  given  tree  the  number  of  that  tree  If 
treenum  values  at  end  nodes  of  a  pivot  eligible  arc  agree,  the  arc  is 
placed  on  LI  Otherwise,  it  is  placed  on  L2  PROCNET  repeats  Step  2 
of  t  tie  primal  simplex  algorithm  for  all  eligible  arcs  from  LI  before 
updating  in  Step  (i  The  length  of  LI  was  set  at  50  and  the  length 

of  L2  was  set  at  30  After  all  pivots  from  LI  have  been  made,  up  to  15 
of  the  best  pivots  from  L2  are  made  following  the  logic  of  | 1 9 ) 

Pure  network  pivots  are  accomplished  following  the  same  procedure 
used  in  a  pure  network  algorithm  For  processing  network  pivots,  yg  in 
(2  10)  is  computed  using  LA05B.  If  l  denotes  a  processing  column  such 
that  [yjjlj  0,  then  PROCNET  flags  basis  trees  containing  processing 


nodes  corresponding  to  column  1.  In  computing  yj  by  means  of  the 
reverse  thread  in  (2  11),  only  trees  which  are  flagged  are  traced 
Since  processing  columns  do  not  leave  the  basis  until  a  reinversion 
occurs,  only  yj  values  are  used  in  the  ratio  test. 

All  parameter  settings  for  PROCNET  mentioned  in  this  section 
remained  fixed  at  these  values  for  the  testing  described  in  Section  5 

8  Computational  Results 

I  n  this  studv  .  test  problems  were  solved  bv  MI’SX  370  and  PROCNET 
Testing  was  done  on  the  IBM  3081-1)  at  the  University  of  Texas  As 
previously  noted.  MPSX'370  is  an  assembly  language  program  while 
PROCNET  is  written  in  FORTRAN  PROCNET  was  compiled  using  the  FORT VS 
compiler  with  optimization  level  3.  and  it  maintains  all  real  values 
using  double  precision  The  execution  times  for  both  codes  are 
reported  in  central  processor  seconds  These  times  do  not  include 

input  or  output  with  the  exception  that  one  line  of  output  (iteration 
log)  is  produced  by  MPSX/370  each  iteration  This  small  amount  of 
output  has  a  negligible  >■  f  f  e<  t  on  the  overall  <  ompa  fisnii  of  t  tic1  two 
codes 

The  CRASH  and  PRIMAL  modules  of  MPSX '370  were  employed  in  solving 
the  test  problems  To  he  comparable  with  PROCNET.  the  reduced  cost 
tolerance  (XTOLDJ)  of  MPSX/370  was  set  to  10  0  It  was  necessary  to 
set  ttie  feasibility  tolerance  ( XTOLV )  to  10  since  the  default  value 
of  10  kept  MPSX/370  from  reaching  feasibility  on  the  test  problems 
All  other  parameters  for  MPSX/370  were  set  to  default  values 


Parameters  used  by  PROCNET .  in  addition  to  those  previously 
discussed,  are  provided  next  The  reduced  cost  and  feasibility 

_  r. 

tolerances  were  both  set  to  10  Big-M  costs  on  artificial  arcs  were 

set  to  99999,  while  pivots  with  minimum  ratio  less  than  10-10  were 
treated  as  degenerate  In  LA05A  and  LA05C ,  pivot  elements  less  than 
0  1  times  tile  largest  element  in  the  pivot  row  were  excluded  Default 
values  were  used  for  other  LAOS  parameters 

The  class  of  allocation  processing  (AP)  network  problems  was  used 
for  testing  These  problems  have  a  dual  block  angular  form  where 
subproblems  corresponding  to  diagonal  blocks  are  transportation  prob¬ 
lems  and  coupling  columns  are  processing  columns  This  class  of  test 
problems  was  also  used  in  [4],  [5] 

Problem  data,  with  the  exception  of  total  supply,  was  randomly 
generated  All  constraints  were  formulated  as  equalities.  As  the 
problems  are  generated,  a  feasible  flow  is  created  The  capacity  of 
each  arc  with  finite  capacity  is  set  to  a  parameter  n  times  the 

feasible  flow  generated  for  that  arc  Total  supply  was  set  at  10000 
fui  all  problems  Two  cost  ranges,  A  and  B.  were  used  Cost  range  A 
has  costs  on  allocation  arcs  in  the  range  100  to  150  and  other  arc 

costs  in  the  range  1  to  100  Cost  range  B  has  allocation  arc  costs  m 
the  range  1  to  100  with  other  arc  costs  ill  the  range  100  to  -1 

Because  a  machine  dependent  ( CDC )  random  number  generator  was 

employed  in  the  previous  studies  ( 4  j .  [5],  we  were  not  able  to  include 

the  problems  from  those  studies  here  The  test  problems  solved  are 

similar  to  those  of  [5],  however,  problems  with  more  processing  columns 
and  problems  with  variable  ^  values  are  included  here 


Test  problem  data  is  given  in  Table  1,  Problems  in  the  groups 


1-4,  5-8,  ...  25-28,  have  the  same  network  topology  Also,  groups  5-8 
and  25-28  differ  only  in  their  cost  ranges.  For  problems  4,8,  ,28, 
the  value  of  ^  was  randomly  selected  from  the  interval  ( 1.1, 2.0]  for 
each  arc  Computational  results  for  these  problems  are  reported  in 
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Table  2.  The  count  of  iterations  for  both  codes  begins  with  the  first 


pivot  after  the  start  (CRASH)  procedure  The  number  of  basic  proces¬ 
sing  columns  at  optimality  is  obtained  by  PROCNET.  This  number  pro- 

T 

vides  an  estimate  of  the  dimension  of  the  matrix  Q  near  optimality. 
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The  ratio  of  total  solution  time  for  MPSX/370  to  that  of  PROCNET  is 
1 1  46  Average  time  per  iteration  for  MPXS/370  is  0.0638  sec/i terat i on 
while  for  PROCNET  it  is  0  0283  sec  i\  terat  l  on  The  ratio  of  these 

average  times  per  iteration  is  about  2  25  The  larger  ratio  for  total 
solution  time  can  be  attributed  to  the  superiority  of  the  PROCNET  start 
procedure  and  pivot  strategy  over  that  of  MPSX '370  on  the  problems 
tested 

The  efficacy  of  the  PROCNET  start  procedure  and  pivot  strategy  is 
strongly  dependent  on  /i  as  shown  in  Table  3  The  ratio  of  total 

MPSX/370  solution  time  to  that  of  PROCNET  decreases  from  17.00  for 
^=1  .  1  down  to  8  45  for  These  results  show  that  PROCNET  is 

especially  effective  on  tightly  capacitated  problems.  Average  pivot 
times  for  the  processing  network  code  are  remarkably  stable  as  /j 
varies  Apparently,  the  tendency  to  a  smaller  number  of  basic 
processing  columns  at  optimality  when  n~\  1  is  offset  by  a  smaller 
percentage  of  pure  network  pivots.  On  the  other  hand,  it  seems  that 
the  larger  percentage  of  pure  network  pivots  for  PROCNET  on  problems 
4,8.  ,28  results  in  a  somewhat  smaller  average  pivot  time 

Table  3  Performance  Values  vs  /c 


Problems  ^  PROCNET  % 

PROCNET 

MPSX/370 

MPSX/370  to 

PROCNET 

Pure  Net¬ 

Avg  Pivot 

Avg  Pivot 

Avg  Pivot 

Ttl  Time 

work  Pivots 

Time 

Time 

Time  Ratio 

Ratio 

1,5,. 

.25  1.1  17.3 

0 . 0290 

0.0572 

1 . 97 

17.00 

2,6, 

.26  2.0  20.3 

0  0295 

0.0662 

2.24 

11.95 

3,7.  . 

.27  «  21.3 

0.0290 

0.0671 

2.31 

8  45 

4.8,. 

. ,28  1 . 1-2.0  28  6 

0 . 0269 

0.0637 

2.37 

11.87 

In  Table  4,  variation  in  code  performance  with  the  number  of  pro¬ 
cessing  columns  p  is  given.  Except  for  p=50 ,  average  pivot  times  for 
PROCNET  increase  with  p  while  the  percentage  of  pure  network  pivots  is 
roughly  decreasing  More  nonzeros  per  processing  column  is  perhaps  the 
major  factor  which  makes  the  problems  with  p=50  difficult  for  MPSX/370. 
We  note  that  for  problems  5-8,  the  MPSX/370  to  PROCNET  total  time  ratio 
is  20.7  The  problems  with  p=100  and  p=150  are  easier  than  problems 
with  larger  p  values  for  both  codes  although  MPSX/370  to  PROCNET  total 
time  ratios  are  down  A  possible  explanation  is  that  the 
transportation  subproblems  have  a  somewhat  different  structure  This 
consists  of  fewer  origins,  more  destinations  and  more  arcs  per  origin 
than  those  with  larger  p  values  Also,  when  p  equals  100,  the  number 
of  problem  rows  is  lower  From  the  results  of  Table  4.  we  conclude 
that  PROCNET  remains  very  efficient  for  problems  with  up  to  250  proces- 
s ing  columns . 


Table  4  Performance  Values  vs  Number  of  Processing  Columns 


'rob  1  ems 

1'rot  ess 

PROCNET 

PROCNET 

MPSX/370 

MPSX  370  t  o 

PROCNET 

C  o 1 umn  s 

Pure  Net¬ 

Avg  Pivot 

Avg  Pivot 

Avg  Pivot 

Tt)  Time 

<p) 

work  Pivots 

Time 

T  ime 

Time  Ratio 

Ratio 

1  4 

10 

42  5 

0  0105 

0  0350 

3  33 

8  38 

-8.25-2 8 

50 

21  6 

0  0278 

0  0727 

2  62 

14  32 

9-12 

1  00 

30  3 

0  0172 

0  0293 

1  70 

5  83 

13-16 

1  50 

31  1 

0  0263 

0  0507 

1  93 

7  51 

17-20 

200 

10  2 

0  0401 

0  0706 

1  76 

13  27 

21-24 

250 

12  0 

0  04  1  1 

0  0648 

1  58 

9  58 

6.  Conclusion 


The  new  version  of  our  pure  processing  network  code.  PROCNET, 
derives  part  of  its  efficiency  from  a  working  basis  of  variable 
dimension  An  extension  of  the  LAQ5  subroutines  has  been  described, 
and  this  extension  is  used  in  PROCNET  to  maintain  the  working  basis 
PROCNET  substantially  outperforms  MPSX/370  in  both  time  per  pivot  and 
total  solution  time  on  problems  with  up  to  250  processing  columns  The 
faster  time  per  pivot  of  PROCNET  indicates  that  it  will  be  useful  as 
part  of  a  mathematical  programming  system  while  the  faster  overall 
solution  time  shows  that  it  can  be  used  to  advantage  in  stand  alone 
fashion  Since  PROCNET  is  an  all-FORTRAN  code,  it  is  highly  portable 
as  well 
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